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We study the properties of the outgoing gravitational wave produced when a non-spinning black 
hole is excited by an ingoing gravitational wave. Simulations using a numerical code for solving 
Einstein's equations allow the study to be extended from the linearized approximation, where the 
system is treated as a perturbed Schwarzschild black hole, to the fully nonlinear regime. Several 
nonlinear features are found which bear importance to the data analysis of gravitational waves. 
When compared to the results obtained in the linearized approximation, we observe large phase 
shifts, a stronger than linear generation of gravitational wave output and considerable generation of 
radiation in polarization states which are not found in the linearized approximation. In terms of a 
spherical harmonic decomposition, the nonlinear properties of the harmonic amplitudes have simple 
scaling properties which offer an economical way to catalog the details of the waves produced in 
such black hole processes. 
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I. INTRODUCTION 

A sensitive array of detectors is searching for the grav- 
itational waves produced in violent astrophysical scenar- 
ios. Detecting and understanding the information car- 
ried by these waves will be crucial to elucidate many 
poorly understood phenomena in astrophysics and cos- 
mology. In particular, the formation of black holes and 
their interaction with the surrounding media might pro- 
duce gravitational waves detectable by space or earth 
based observatories, such as LISA and LIGO (see for in- 
stance 0, 0i 13)- It is therefore important to consider 
these kinds of scenarios and produce reliable waveform 
estimates which can be used as input for data analysis 
techniques. 

In the case of LISA, its frequency band will be sensitive 
to systems involving supermassive black holes. While 
the formation process of these is not well understood, 
different working models suggest the following mecha- 
nisms: formation by a series of smaller black hole merg- 
ers; growth from smaller holes by accretion or collapse of 
massive gas accumulations Q. Irrespective of the mech- 
anism, after some stage it would lead to the scenario of a 
massive black hole which is being strongly perturbed. For 
earth based detectors, systems of interest include merg- 
ers of stellar mass black holes and/or neutron stars. The 
early stages of these systems can be fairly accurately de- 
scribed by post-Newtonian approximations (e.g. 
while the late stage can be described by linearized pertur- 
bations off a single black hole spacetime. The intervening 
strong field stage can only be handled accurately via nu- 
merical simulations of Einstein equations which, to date, 
are not yet available (see 0, @ for a review on the status 
of efforts in this direction). 

Irrespective of whether the massive black hole formed 
from accretion of matter or from mergers of smaller black 



holes, quite useful information can be obtained by ex- 
ploiting the fact that at some stage the system can be 
treated as a perturbed, single black hole spacetime (as 
was demonstrated in |9lll0l|). 

Therefore, a numerical code that can stably deal with 
generic single black hole spacetimes can be useful in the 
study of nonlinear disturbances, probing whether there 
are robust features which are reflected in the waveforms 
produced. If these features are present, the information 
can be incorporated in data analysis (see for instance 
)- which would be of considerable value until more 
refined waveforms can be obtained. 

A mature characteristic evolution code, the PITT null 
code, is capable of handling generic single black hole 
spacetimes [Lj . The code computes the Bondi news func- 
tion describing radiation at future null infinity X + , which 
is represented as a finite boundary on a compactified nu- 
merical grid. The Bondi news function is an invariantly 
defined complex radiation amplitude N = N§ + iN®, 
whose real and imaginary parts correspond to the time 
derivatives and dth® of the "plus" and "cross" po- 
larization modes of the strain incident on a gravitational 
wave antenna. An alternative approach to calculating ra- 
diation in black hole spacetimes might be based upon the 
Cauchy formulation, which has recently undergone signif- 
icantly improvement p^ . ll5lllrill7| . In the present work, 
we employ the characteristic code, developed in |13| and 
further refined in [Isj . to study the nonlinear response 
to the scattering of gravitational radiation off a non- 
spinning black hole (the case for spinning black holes 
being deferred for future work). Although the gravita- 
tional waves will be weak when they reach a detector, 
they are produced in a region where nonlinear effects are 
important. 

In this work, we investigate in a simple setting how 
such nonlinearities affect the produced waveforms. Con- 
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siderably more work certainly lies ahead but already im- 
portant effects are noticed. Namely, by perturbing the 
black hole with a pulse with a single (£, to) mode of am- 
plitude A we observe: 

• There is generation of additional modes whose am- 
plitudes scale as precise powers of A. This suggests 
an economical way of producing a waveform catalog 
based upon combining the linearized results with a 
modest number of (non-linear) simulations. 

• There is significant phase shifting and frequency 
modulation in the nonlinear waveforms. 

• There is nonlinear amplification of the gravitational 
wave output. 

• There is a significant generation of radiation in po- 
larization states not present in the linearized ap- 
proximation. 

Throughout this work to illustrate the above men- 
tioned effects we standardize the input pulse to (£ = 2, 
to = 0) and (£ — 2, to = 2) quadrupole modes, although 
the simulations can be carried out with any input data. 
Our main objective here is to demonstrate the effective- 
ness of the characteristic code in providing a complete 
analysis of the mode coupling in the outgoing waveform. 
The robustness observed in this simple case opens the 
door to performing in-depth analysis of more general in- 
put pulses. For example, a localized ingoing pulse of 
gravitational energy, patterned after the distortion in the 
gravitational field around a localized distribution of mat- 
ter, could be used to mimic the effects of the infall of 
matter onto a black hole until more realistic calculations 
are achievable. 

Mode mixing and the transition from the linear to the 
nonlinear regime have previously been studied by Allen et 
al. 0| and Baker et al. (2*3 m a Cauchy evolution frame- 
work. They extract nonlinear waveforms by matching the 
perturbative and nonlinear solutions on a worldtube at 
r = 15M, where M is the unperturbed black hole mass. 
Our work is complementary to theirs in the sense that 
we use a characteristic formalism to carry out the evolu- 
tions. A characteristic approach offers greater flexibility 
and control in prescribing initial data. In the Cauchy 
approach, elliptic constraints must be solved in order to 
provide initial data. In order to simplify the constraint 
problem, Allen et al. [l9| use time-symmetric initial data, 
which intrinsically contain equal amounts of ingoing and 
outgoing radiation. In our characteristic approach, we 
prescribe initial data on a pair of intersecting null hyper- 
surfaces, one of which is ingoing and the other outgoing. 
Data on each null hypersurface can be prescribed freely 
to represent a pulse with arbitrary waveform. Further- 
more, the data on an outgoing null hypersurface can be 
identified with an ingoing wave; and data on an ingoing 
null hypersurface, with an outgoing wave. There is no 
comparable way to prescribe purely ingoing or outgoing 
Cauchy data. 




FIG. 1: The physical setup for the scattering problem. A star 
of mass M has undergone spherically symmetric collapse to 
form a black hole. The ingoing null worldtube M lies outside 
the collapsing matter. The metric inside M (but outside the 
matter) is vacuum Schwarzschild. Outside of A/", data for an 
ingoing pulse is specified on the initial outgoing null hyper- 
surface J~ . As the pulse propagates to the black hole event 
horizon TC + , part of its energy is scattered to X + . 



A characteristic treatment of mode mixing in the ax- 
isymmetric case (to = 0) has previously been undertaken 
by Papadopoulos |2lj | , using the initial axisymmetric ver- 
sion of the PITT code. Papadopoulos carried out an illu- 
minating study of the propagation of an outgoing pulse 
by evolving it along a family of ingoing null hypersur- 
faces with outer boundary at r = 60M. The evolution 
is stopped before the pulse hits the outer boundary in 
order to avoid spurious reflection effects and the radia- 
tion is inferred from data at r — 20M. As in the Cauchy 
treatment of Ref 's 0, [2"3 , gauge ambiguities arise from 
reading off the radiation waveform at a finite radius. 

In the present work, we study the scattering of both 
axisymmetric and non-axisymmctric ingoing pulses by 
evolving along outgoing null hypersurfaces. The physical 
setup is described in Fig. ^ The outgoing null hypersur- 
faces extend to future null infinity X + on a compactificd 
numerical grid. Consequently, there is no need for cither 
an artificial outer boundary condition or an interior ex- 
traction worldtube. The outgoing radiation is computed 
in the coordinates of an observer in an inertial frame at 
infinity, thus avoiding any gauge ambiguity in the wave- 
form. Although the technical setup is very different from 
the previous treatments figL f20L l21| . the simulations dis- 
play several qualitatively similar features (see Sec. lVI|) . 

In Sec. m we briefly summarize the characteristic for- 
malism on which the paper is based. In Sec. IIIII we give 
the geometrical setup for initializing and setting bound- 
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ary conditions for the simulations. In Sec. IIVI we cali- 
brate the accuracy of the full nonlinear code against re- 
sults obtained with independent linearized characteristic 
codes |2l[23| that solve the Teukolsky equations plE^ 
for perturbations about Schwarzschild spacetime. We 
show, in the regime where linear perturbation theory is 
valid, that the full nonlinear code reproduces the Bondi 
news function obtained from the perturbative codes to 
second order accuracy. 

Simulations of nonlinear mode-mode coupling are de- 
scribed in Sec.0 The nonlinear effects on the waveform 
of the outgoing radiation can be cleanly separated out 
by comparison with the waveform produced by the lin- 
earized code. The simulations reveal several nonlinear 
features of the waveform that have potential importance 
to the design of templates for data analysis. These are 
summarized in Sec. I VII 



II. THE NONLINEAR SCATTERING PROBLEM 

Following 0, 0, |2(| , we use coordinates based 
upon a family of outgoing null hypersurfaces, letting u 
label these hypersurfaces, x A (A = 2,3), label the null 
rays, and r be a surface area coordinate, such that in the 
x a = (u, r, x A ) coordinates the metric takes the Bondi- 
Sachs form HI 



spherical harmonics, is given in App. [S] For details, see 
Ref. [13. 

The vacuum Einstein equations form (i) a hierarchy 
of null hypersurface equations for the radial derivatives 
of P, U and W (and also for auxiliary variables v : k, B 
and Q used to reduce the equations to first order differ- 
ential form in the angular coordinates), (ii) a complex 
evolution equation for the metric function J and (iii) a 
set of subsidiary equations which are satisfied by virtue 
of the other equations if they are satisfied on a null or 
timelike worldtube. In our case, a null worldtube forms 
the inner boundary of the computational domain, and we 
solve the equations induced there in the form discussed in 
Ref. • Most of the discussion in the present paper will 
deal with an "ingoing" pulse, in which case Schwarzschild 
data are posed on the inner boundary, the subsidiary 
equations are then solved automatically. The explicit 
form of the equations (i) and (ii) is given in Eq's. (16)- 
(26) of Ref. pjj. 

Given initial data J on an outgoing null hypersurface 
extending to X + and boundary data on an inner timclike 
or null world tube, the PITT null code carries out the 
future evolution of the spacetime. The evolution extends 
to I + where the Bondi news function N is calculated. In 
the present work, the computation of the news is carried 
out with enhanced accuracy by the new method described 
in App. |B] 



ds 2 



' \ -) r 2 h AB U A U B 



r 

2r z h AB U B dudx A + r z h A Bdx A dx 



du 2 - 2e 2p dudr 
(2.1) 



Here W is related to the Bondi-Sachs variable V by 
V = r + W and det(h A B) — det(q AB ), with q AB a stan- 
dard unit sphere metric. We represent q A B in terms of a 
complex dyad q A satisfying q A q A = 0, q A q A = 2, q A — 
q AB q B , with q AB q B c = & A and q AB = \ {q A q B + qAqB)- 
The angular coordinates x A — (q,p) used in the code are 
based upon a complex stereographic coordinate z = q+ip 
for the unit sphere metric 



q AB dx dx 



1 + q 2 + p- 



idq 2 + dp 2 ), 



tanfe^ (see 



with two patches used to cover the sphere. In the north 
patch, the stereographic coordinate is related to stan- 
dard spherical coordinates (0, <f>) by zn = t£,n - pl ^ 
App. . 

We represent tensors on the sphere by spin-weighted 
variables |30l |. Thus the conformal metric h AB is rep- 
resented by the complex function J = \h AB q A q B , and 
the real function K = \h AB q A q B , where K 2 = 1 + J J 
on account of the condition det(h AB ) = det(q AB ). The 
metric functions U A are similarly encoded in the com- 
plex function U = q A U A 0,H||. Angular derivatives 
are expressed in terms of S and 5 operators. A brief de- 
scription, along with our conventions for spin-weighted 



III. DATA FOR THE SCATTERING PROBLEM 

The nonlinear simulations described in Sec. [V] model 
the scattering of an ingoing pulse in the region exterior to 
a spherically symmetric collapsing star of mass M, as de- 
picted in Fig. ^ (the mass M will always be scaled to 1 in 
our simulations). The free initial data for the pulse con- 
sist of the metric function J on the initial outgoing null 
hypersurface, denoted by J~ . In addition, Schwarzschild 
null data is given on the interior null hypersurface Af, 
which is causally unaffected by the pulse. The evolution 
code then provides the news function at 2 + . For compu- 
tational simplicity, we simulate a mathematically equiv- 
alent problem by considering the Kruskal continuation 
of the Schwarzschild interior so that the inner bound- 
ary is extended to the ingoing branch of the r = 2M 
null hypersurface, which we denote by Tl~ . In our met- 
ric ansatz ()2.1|l Schwarzschild data on ~Hr correspond to 
setting /3 = 0, U = 0, W = —2M. A schematic view 
of this setup, in the linear regime, is provided in Fig. [5] 
(Fig-E]shows the compactified exterior Kruskal quadrant 
which includes T~ as well as T + ). In the nonlinear simu- 
lations, the data on Ti~ is set to the Schwarzschild data 
J = 0, which implies that the intended Schwarzschild 
data is also induced on N. In the code calibration tests 
considered in Sec. IIVI we also consider outgoing pulses 
generated by data on H~ . 

We used two different types of evolution coordinates, 
u and u. For the calibration runs we chose as our evo- 
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lution coordinate u which is an afhne parameter along 
the null generators of 7i~ . For the remaining runs we 
chose as our evolution coordinate u which is the stan- 
dard Schwarzschild retarded time on H.~ (u = t — r*). 
Later, when analyzing gravitational wave signals, we will 
use Bondi time - an afhne parameter u along the null gen- 
erators of X + . For Schwarzschild, Bondi time coincides 
with the standard Schwarzschild retarded time u = u, 
and u = -Mexp(-u/AM). 

We prescribe the input data for the pulse in the form 
J = f{x)2Rim, where x — r for initial data on J~ 
or x = u for data on TL~ and s Rim are spin-weighted 
spherical harmonics with real potentials (see App. lAjl . 
The s Re.m are linear combinations of the standard s Yi m 
and form a complete basis for smooth spin- weight s func- 
tions. In linear theory a particular 2 Rim mode in the 
data produces a Bondi news function that contains only 
that mode. Conversely, data constructed from the stan- 
dard spin- weighted spherical harmonics 2^im would pro- 
duce both the 2Ye m and 2 Ye - m modes in the Bondi news. 
An explanation for this linear mode-coupling and explicit 
formulae for the s Rim are given in Add. lAl 

The function / determines the radial profile of the 
pulse. As for its initial angular dependence, we restrict 
our attention here to input data with I = 2 and m = or 
772 = 2. In the regime where the linear approximation is 
valid, the Bondi news function N will also be of this form 
but nonlinear effects lead to the presence of higher order 
modes. Parity and reflection symmetry implies that only 
even £ and even m > modes be generated by such non- 
linear effects. 

The input pulse is specified in terms of the function 



m = a 



2 {x Xjnin) {x Tl 



(x r , 



v 
1 1 



(3.1) 



where A controls the amplitude and 77 controls the steep- 
ness. For outgoing pulses, we prescribe boundary data 
on with the specific profile 



J\ H -(u,x A ) = f(u) 2 R 2 , 



(3.2) 



with 77 = 6, for u m i n < u < u max ; elsewhere we set 
J = 0. For an ingoing pulse, we give data on the initial 
outgoing null hypersurface of the form 



J(r,X A ) = f( X ) 2 R2m 



(3.3) 



for x m in < x < x max , where x — r/(2M + r), with cither 
72 = 3 or 77 = 6; elsewhere we set J = 0. We also set 
J = on the inner boundary 7i~ . 



IV. CALIBRATION OF THE NONLINEAR 
CODE AGAINST PERTURB ATIVE SOLUTIONS 

We calibrate the Bondi news function obtained from 
the PITT nonlinear code against values obtained from 
a characteristic perturbative code [22T [2j| based upon 



the Teukolsky equations for the Newman-Penrose quanti- 
ties [3^ 4>o = C a bcdJ a, m b l c 'm d and $4 = C a b c d'n a fh h h c fh d . 

We fix the null tetrad by setting l a = — V a tt and 
772 a = q a /\/2r, where u = t — r* in the background 
Schwarzschild spacetime [r* = r + 2Mlog(r/2M — 
1) is the Regge- Wheeler tortoise coordinate]. Rather 
than evolving the spin- weighted quantities ij)Q and ^4, 
we evolve equivalent spin-zero potentials (see App. 
rescaled by appropriate factors of r. These evolution 
variables are F4, where S 2 ^ = rip4, and Fq, where 
B 2 F = (1 - 2M/r) 2 r 5 4. Both F 4 and F are con- 
structed to be hnite on I + . 

The perturbative variable F± has a complicated depen- 
dence on the Bondi metric variables. However, on H~ , 
where the data for the outgoing pulse is prescribed, the 
expression for F4 reduces to 

„,2 



where u = — Mexp(— u/AM) is an afhne parameter along 
the null generators of "Hr and j is the spin-zero potential 
for J, i.e. J — <5 2 j. On an outgoing null hypersurface, 
where the data for an ingoing pulse is prescribed, the 
perturbative variable Fq has the simple dependence on j 
given by 



Fo = ^r(r-2M) 2 d r (r 2 d r j). 



(4.2) 



For input data consisting of an outgoing pulse, where we 
specify compact support boundary data for J on HT and 
set J = on J~ , we directly obtain equivalent boundary 
data for F4 via Eq. I|4.1|l . For input consisting of an in- 
going pulse, where we specify the compact support data 
on J~ , we directly obtain the equivalent boundary data 
for Fq via Eq. (|4.2|l . For the calibration tests presented 
here we choose input data with (I = 2, 772 = 0) spher- 
ical harmonic dependence for which F4 — /4(u, r)i? 2 
[Rem = oRim] and F = fo(u,r)R 2 o- The perturbative 
code provides solutions for f 4 and fn- Details of the evo- 
lution code for f 4 can be found in |22(; and details for /o, 
in H|. 

In the perturbative approach, the Bondi news function 
is either obtained directly from F4 on I + or indirectly 
from Fq The Bondi news function is given by 



N(u) 



F t 



4 1+ 



(t)dt, 



(4.3) 



where the integral in Eq. (|4.3|l is over Bondi time. In the 
simulation of an outgoing pulse, the input data for F4 on 
"Hr is set to zero until the start of the evolution at u = 
Umin- In this case we calculate the Bondi news through 
a second-order accurate midpoint-rule integration of F4. 

In the simulation of the scattering of an ingoing pulse, 
the computation of the news function is more compli- 
cated. For the present case of an (£ — 2, 777 = 0) lin- 
earized mode, the news function is related to Fq by 



N - 



M 



1 



duN = -Q 2 diF \ I+ 




(4.4) 
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Care must be taken in initiating the integration of 
Eq. (|4.4I) because the Bondi news function N does not 
vanish at the initial evolution time. We need an accu- 
rate initial value for N, which we obtain by integrating 
Eq. (|4.2() to obtain the metric variable J. In the linear 
regime, the Bondi news function can be calculated from 
the value of J and its radial derivative at 1 + ^| • I n the 
present case of an (£ = 2, m — 0) linearized mode, 

N = -^(3J + r 2 d r duJ)\ x + (4.5) 

(a generic formula for the linearized News as a function 
of J|j is given in Eq. H5.2fl ). We use Eq. I|4.5|l to obtain 
the news function on the initial slice and also to compute 
it at later times as a check that it produces the same 
values as the integration of l|4.4|l . 

Because the perturbative code is one-dimensional, 
computations on very large radial grids can be carried 
out to provide effectively exact solutions for calibrating 
the error in the Bondi news function calculated using the 
nonlinear code. Such a check must be carried out in the 
range of validity of the linear approximation, and the er- 
ror in the perturbative calculation must be sufficiently 
smaller than the discretization error of the nonlinear null 
code. For this purpose, simulations with the perturba- 
tive code were carried out with 16001 radial grid points 
for the outgoing pulse. For the ingoing pulses (for which 
the perturbative code is less accurate), we used 2501, 
5001, and 10001 radial grid points, and then performed 
a Richardson extrapolation to obtain high accuracy. 

The initial and boundary data given by Eq's. I|3.2|) and 
(|3.3|) are compact pulses with relatively steep profiles, 
which require comparatively large grid sizes to attain 
good resolution. The specific form of the profiles lead to 
C 3 differentiability of the associated perturbative quan- 
tities i*4 and Fq, a requirement which ensures that the 
news function can be obtained from the perturbative cal- 
culation to second order accuracy. In all tests we scale 
the initial mass of the system to M = 1. 

A. Propagation of outgoing pulses 

In our first battery of tests, which provide stringent 
tests of the ability of the code to carry radiation away 
from the horizon Ti.~ to null infinity X + , we prescribe 
an arbitrary perturbation on the horizon while setting 
the perturbation of the initial outgoing null cone to zero. 
The horizon data for J are given by Eq. I|3.2[l with 
Umin = —1.99, u rnax — —1.5, and A — 2.6 ■ 10~ 5 . This 
corresponds to a relatively short pulse of Au ss 1.13M in 
the background Schwarzschild retarded time. The corre- 
sponding data for F4 are obtained by applying Eq. (|4.1|) 
to Eq. i)3. 2f> . We measure the error in the Bondi news 
function N code with the L m norm \\N code - N per t\\oo re- 
stricted to the interior of the equator on the stereographic 
patches. Our primary concern is to check convergence 
and the evolution was stopped at u — —1.5 when the 
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FIG. 2: The rescaled norms £ = (i/6) 2 \\N code — N pert \\oo 
(versus time) of the Bondi news for the compact support out- 
going pulse. 



horizon data vanishes. The perturbative news function 
was obtained using 16001 radial grid points and an initial 
time-step of du = 8 • 10~ 5 . 

The nonlinear characteristic code uses a uniform ra- 
dial grid of n x points on the compactified coordinate 
x = r/(R + r), where R = 2M is the initial horizon ra- 
dius. The coordinate range isl/2<x<l, with x = 1/2 
at Hr and % = 1 at I + . We introduce two additional 
ghost zones inside the horizon in accord with the start-up 
algorithm described in |27| . The angular grid is also uni- 
form, with n q grid points spanning the range (— qs,qs) 
in both the q and p directions. We set qs = 1.2 in order 
to provide a finite overlap between the two stereographic 
patches. In the convergence tests we take n q = 12t + 5, 
n x = 30z + 3, with i taking integer values from 6 to 11. 
We keep the time step du constant during each run, at 
du = 5.6 • 10~ 3 /i, below the smallest value required for 
the CFL condition to be satisfied and scaling with i to 
guarantee convergence. 

Since the code is second order accurate, and the lin- 
ear approximation holds, the error norm should be in- 
versely proportional to the square of the grid size, i.e. to 
i 2 . Figure |21 shows the rescaled norm of the error 
for the above 6 grids. Note the perfect overlap at early 
times (prior to u — —1.65). The later errors scale ap- 
proximately with i 2 but show definite deviations. Runs 
with smaller amplitude indicate that this error is not a 
nonlinear effect. The error appears to be due to accu- 
mulation of higher order truncation error, mixed with 
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FIG. 3: The news and norm £ = || JVcode — N per t\\oo 
(rescaled by a factor of 2000) for the compact support outgo- 
ing pulse. 



FIG. 4: The rescaled L 2 norm £ = (i/6) 2 1| A^ code - N pert \\ 2 
versus grid-parameter i for the ingoing pulse test. Second 
order convergence is confirmed by the reasonable overlap of 
the norms. 



some smaller roundoff error. Nevertheless, the calcula- 
tion is still extremely accurate including up to the end 
of the run, when deviation from strict second order con- 
vergence is more noticeable. This is evident from Fig. |3J 
which shows the news function and the error norm 
(multiplied by 2000) for the finest grid. The error is 
at least 200 times smaller than the news function; thus 
controlled errors of less than one-half of a percent are 
achievable at the larger resolutions. 



B. Scattering of ingoing pulses 

These tests correspond to the realistic situation of the 
scattering of gravitational radiation which has been cre- 
ated in the near field of a black hole. We place a com- 
pact support pulse on the initial null hypersurface and 
set the boundary data on H.~ to zero. The initial data 
were given by Eq. (|3.3|l . The setup is described in Fig.^ 
The perturbation calculation was carried out using the 
evolution algorithm for Fq described in |23|. For this 
test we choose x m i n = 0.6, x max = 0.8, n = 6, and an 
amplitude A of 2.6 • 10~ 6 . The corresponding pulse ini- 
tially extends from r = 3M to r = 8M. The nonlinear 
runs were performed with the same grid parameters as 
in Sec. lIV A1 with the grid sizes again determined by the 
integer i ranging from 6 to 11. The time step was set to 
du = —ub ■ 10~ 4 /i, where the factor of u ensures that 
the CFL condition remains satisfied (by approximately 



a factor of 1 /4) . Figure 01 shows the rescaled Li norm 
of the error £ = (i/6) 2 | N co d e — N pert \\2 for the above 6 
grids. The excellent overlap of the norms confirms that 
the nonlinear null code is second order convergent. 

The second order convergence of the news from the 
ingoing pulse begins to break down at u = 15M as a 
consequence of our choice of coordinate system. Fig. 
shows a plot of Schwarzschild background spacetime with 
the timelike curves corresponding to the worldlines of 
the first few radial gridpoint (the ingoing null curve TL~ 
is the innermost gridpoint). An initially compact pulse, 
bounded by the ingoing null lines K\ and K2, is placed on 
J~ . At the time level J + the ingoing pulse will occupy 
a region containing only a few radial gridpoints. Con- 
sequently, a well resolved pulse on the initial slice will 
cease to be resolved in a finite amount of time. Radial 
resolution of the ingoing pulse begins to break down at 
u = 15 M. This breakdown in resolution happens rel- 
atively quickly because we chose a grid that is uniform 
in the x coordinate. The evolutions could be extended 
by using adaptive mesh refinement techniques recently 
developed for characteristic codes [3i| , and by different 
choices for the compact radial coordinate. (Note that the 
loss of resolution happens in the radial direction, not in 
the angular ones, so that the techniques in j33| could be 
readily applied to this case. We will revisit this point in 
Sec. ED 
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FIG. 5: The ingoing pulse on a Schwarzschild background. 
The timelike curves (along with TL~) are the worldlines of 
the first 5 radial gridpoints (the first gridpoint is located at 
TC~). The r = const lines begin at past timelike infinity i~ 
and terminate at future timelike infinity i + . An ingoing pulse 
well resolved on J~ is compressed (in the x coordinate) into 
3 gridpoints by the time it reaches J + 



V. ANALYSIS OF MODE-MODE COUPLING IN 
GRAVITATIONAL RADIATION SCATTERING 

The (u, x, x A ) coordinates, which define the computa- 
tional grid of the PITT code [l3L ITsI ] are adapted to the 
geometry of the inner boundary. A gravitational wave de- 
tector would measure the radiation as seen in a distant 
inertial coordinate system, located essentially at I + . Wc 
denote by (uB,y A ) the corresponding inertial Bondi co- 
ordinates on 1 + (where x = 1). Thus for the purposes of 
gravitational wave data analysis the news function must 
be expressed in the form N(ub,U A )- I n the course of 
this work, modules have been added to the PITT code 
to carry out this transformation, which is essential in or- 
der to provide correct time profiles and harmonic mode 
analysis. In the linear regime, (u,x A ) ss (ub,V A ) so that 
these corrections to the news function are of second order 
and can be ignored. That is no longer the case when we 
consider the nonlinearities manifested by mode coupling. 

The "inertial" news function N(ub, y A ) is obtained by 
performing a fourth order accurate interpolation between 
the (it, x A ) and (mb, y A ) grids. It is then decomposed into 
spin- weight 2 spherical harmonic amplitudes Ng m via 
second order accurate integration over the sphere with 
solid angle fl = 4ir, 

N im = I N 2 Rt m dtt. (5.1) 



See App. [S] for further details concerning the spin- 
weighted harmonic decomposition and how the integra- 
tion is carried out on the stereographic patches. 

In linear theory, in a gauge where f3 — 0(A 2 ) and 
U = 0(^4), the Ng m coefficients are related to the Jg m 
coefficients of the metric function J (obtained by apply- 
ing Eq. (|5.1|l to J rather than N) by 

N lm = --r 2 d r d u J em - K ± ' ^{Jt m ). (5.2) 

In linear theory there is no mode coupling in J, and initial 
data containing a single mode produces a Bondi 

news function containing only that mode. 

A. Axisymmetric mode coupling 

We first study mode coupling in the axisymmetric case 
by prescribing initial data as an ingoing [t. = 2, m = 0) 
pulse, as in Eq. (|3.3(l with n — 3, x m { n = .56, and 
Xmax = -8. This choice of initial pulse, extending from 
r = 2.55M to r = 8M, is made for computational econ- 
omy because the nonlinear effects arc weak for r 3> 3M. 
We vary the amplitude from A = 2.6 • 1CP 6 , well in the 
linear regime, to A = 2.6 • 10 _1 , where nonlinear effects 
can be clearly observed. We display results obtained with 
a resolution of 77 x 77 angular grid points and 181 ra- 
dial grid points. This grid size (the size of the base grid 
used in the convergence tests) allows a broad search for 
interesting qualitative behavior with reasonable compu- 
tational time. The news function is computed in inertial 
Bondi coordinates and then decomposed into spin- weight 
2 spherical harmonics. The reflection symmetry and ax- 
isymmetry of the initial data are preserved by the numer- 
ical evolution in the nonlinear regime, so that the only 
non- vanishing amplitudes N4 m °f the scattered radiation 
have even I and m = 0. In addition, these symmetries 
imply that the Ngo are real, which in our conventions 
means that only the © polarization mode is excited. 

The following figures show the time dependence of the 
lowest order mode amplitudes, N20, N40, and Nqq, which 
are excited in the news function by a given input ampli- 
tude A. Figure shows the dependence of N20/A on 
A. From the figure it is clear that nonlinear effects do 
not make N20 deviate significantly from its perturbative 
waveform for A < 2.6 • 10~ 2 . However, at higher am- 
plitudes the effect of nonlinearity is to amplify the wave, 
which could enhance the prospects of detection above the 
level predicted by linearized theory. This nonlinear am- 
plification is primarily quadratic in A. In addition to the 
nonlinear amplification there is also a phase shift which 
can be seen more clearly in Fig. 1111 

Nonlinear effects are more easily seen in the modes that 
vanish in the linearized approximation. Figure [7] shows 
that A/40 scales as A 2 at sufficiently high amplitudes. For 
smaller amplitudes A this quadratic scaling is masked by 
the truncation error introduced in the computation of the 
spherical harmonic decomposition of the news function at 
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FIG. 6: The rescaled coefficient N 20 /AforA = 2.6- 10"\ 2.6- 
10" 1 - 2 ,2.6-l(T 1 ' 5 ,2.6-l(nV--,2.6-10~ 6 . For A < 2.6- 1(T 2 , 
N 2 1 A has negligible dependence on A and the curves overlap. 



the given grid size and, as a result, A40 scales roughly 
linearly with A. 

Figure |H1 shows N^o scales as A 3 at the highest ampli- 
tudes. As might be expected for a higher £ mode, the 
masking effects of truncation error are now more severe 
at the lower amplitudes. The breakdown of scaling be- 
havior in Fig's. El - El results from inaccuracy at times 
beyond u = 20 due to lack of resolution of short scale 
features using the current grid. These short scale fea- 
tures arise inside r = 3M where quasinormal ringing is 
known to originate. 

There are some understandable aspects in these mode 
coupling results. Apart from numerical truncation error 
effects at low amplitudes, we find that the corresponding 
amplitudes A^o scale as A 1 -/ 2 . This is consistent with the 
property of ordinary spherical harmonics 



{Y lm y 



Y, 



<Jt) ■< 



terms with smaller I. 



Thus we expect that an £ — 2j mode arises from order 
j (and higher) nonlinear terms, and hence will scale as 
A? + 0(A : ' +1 ). The theory of how the Einstein equations 
couple these spin-weighted spherical harmonics has not 
been worked out. A full analysis of the possible modes 
produced at a given order would require a computation 
of the appropriate Clebsch-Gordon coefficients for spin- 
weighted spherical harmonics. Finding these coefficients 
is complicated by the fact that there are many ways to 
combine £ and m modes of various spin-weights to pro- 
duce a spin 2 function (e.g. ±Rg m %Rt' m '; \Ri m iRt'm>; 
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FIG. 7: The rescaled coefficient A/40/Afor A = 2.6- 10~ 3 , 2.6- 
10 -4 , • ■ • , 2.6 • 10 -6 , and the rescaled coefficient N40/A 2 for 
A = 2.6 • 10 _1 ,2.6 • 10 _1 - 2 ,2.6 • 10- 15 ,2.6 ■ 10~ 2 . Note that 
at early times N40/A is independent of A for A < 2.6 ■ 1CF 4 . 
For small A, the computed value of A^o is dominated by 
truncation error and is thus proportional to A. For larger 
A, where its computed value is of physical relevance, A^o 
is proportional to A 2 . When A = 2.6 • 10 -3 , the coefficient 
contains a mixture of order O(A) truncation error and order 
0(A 2 ) non linear terms. 



oRim iRi 1 m't ' ' ')■ Rather than work out these coeffi- 
cients we use a more naive approach of looking at the 
possible modes produced by combinations of spin- weight 
harmonics. From this approach we find that the pos- 
sible modes produced by combining an [£ = 2, m = 0) 
mode with itself (i.e. quadratic coupling) are {£ = 4, 
to = 0) and {£ = 2, m = 0) [£ = is not allowed for 
spin 2 fields]. Thus the (£ — 2, to = 0) mode in the news 
would contain linear and higher order terms, whereas the 
{£ = 4, to = 0) mode would contain quadratic and higher 
terms. The above numerical results are consistent with 
these naive expectations. In addition, there is a similar 
behavior in the frequencies of the various £ modes. Fig- 
ure |3| shows the Fourier transforms of the A/20, A/40, an d 
Aeo coefficients for the A = 2.6 • 1CT 1 run. The (£ = 2, 
to = 0) mode shows a strong peak with maximum at 
iv .36; the (£ = 4, m = 0) mode has peaks at uo » .73 
and lu w 0; and the {I = 6, m = 0) mode has peaks at 
u ~ 1.10 and lu « .30. This entire behavior is exactly 
what would arise from the nonlinear power law response 
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FIG. 8: The rescaled coefficient N 60 /A for A = 2.6 
10 -2 , ■ ■ • , 2.6 ■ 10 -6 , and the rescaled coefficient Nqq 
A = 2.6 • 1CT\2.6 ■ 10 _12 ,2.6 ■ 10~ 15 . Note that N aQ /A is 
independent of A for A < 2.6 ■ 10~ 2 . For small A, the value 
of the A^6 o is dominated by truncation error and is thus pro- 
portional to A. For larger A, the coefficient is proportional 
to A 3 . However, the nonlinear terms in Neo are smaller than 
truncation error for A as large as 2.6 ■ 10 -2 . 



/ f / 2 to an I = 2 mode /; i.e. quadratic terms produced 
from / ~ sinwoi contains frequencies 2wo and 0. Simi- 
larly, cubic terms produce frequencies wq and 3wo- 

The dominant contribution to the frequency of the 
[i = 2, m = 0) mode is expected to be the lowest £ = 2 
quasinormal frequency. For an M = 1 black hole the 
lowest I — 2 quasinormal frequency is u> = .373672 [34[. 
However, the A = 2.6 - 10 _1 amplitude pulse should con- 
tribute significant mass to the system and a lower fre- 
quency should be expected, in agreement with the mea- 
sured peak at w w .36. To measure the accuracy of the 
frequencies obtained for these systems we performed a 
similar analysis with the A = 2.6 • 10~ 6 run. In this low 
amplitude limit, the Fourier transform of -/V20 showed a 
maximum at u> ss .37, very close to the frequency of the 
lowest quasinormal mode. 

Of direct importance for designing templates for wave 
detection is the waveform obtained by the net superpo- 
sition of these modes. In Fig.^]we plot the time depen- 
dence of the inertial news function, as measured by an 
observer at the equator in a coordinate system adapted 
to X + , as a function of input amplitude. The figure shows 




FIG. 9: Fourier transforms (after rescaling to similar ampli- 
tudes to facilitate comparison of peak locations) of the N20, 
-/V40, and Neo time profiles for the A — 2.6 • 10 _1 run. The 
(I = 2, m = 0) mode shows a strong peak at to ~ .36, the 
(I — 4, m — 0) mode shows a strong peak at ui ~ .73 and 
lj ~ 0, and the (I = 6, m = 0) mode shows a strong peak at 
lo w 1.10 and u w .30. 



nonlinear amplification of the news function with increas- 
ing input amplitude. The peaks of the waveform differ 
by up to a factor of 2 from the values which would be 
obtained by a linearized calculation. 

In addition there are phase shifts in the location of the 
peaks. These are apparent in Fig. 1111 which compares 
N/A at the equator for the nonlinear case A — 2.6 • 10 _1 
and the linear case A — 2.6T0~ 6 . The difference between 
these two waveforms corresponds to the error that would 
be made by using a linearized calculation to obtain the 
A = 2.6 ■ 10 _1 waveform. In the first oscillation the A = 
2.6- 10 _1 run has the larger wavelength, while afterward 
its wavelength is shorter. 

Accurate knowledge of the phase is very important to 
data analysis for extracting the signal [3{|. When using 
matched filtering over a number of cycles of the wave- 
form, the total integrated error in the phase must be no 
greater than 10% of a cycle. For example, a phase shift of 
2% per cycle from the value programmed into the detec- 
tion template would render the phase information useless 
in 5 cycles. The trend displayed in Fig's. [Trjlandllll is a 
drift in phase of about 15% from the first maximum to 
the first minimum. The nonlinear news oscillates quicker 
in this first half cycle. This trend is reversed in the sec- 
ond half of the first cycle and the net phase shift relative 
to the linear news after the first complete cycle is 2%. 
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FIG. 10: The rescaled Bondi news function N/A on the equa- 
tor. The plot is very similar to Fig. |S| (the plots differ by an 
overall factor which results from the normalization of 2^20 
and phase shifting in the larger amplitude runs) due to the 
dominance of the (t = 2, m = 0) mode. The news is axisym- 
metric and contains only the © mode. 

Beyond the first cycle the nonlinear news exhibits a con- 
sistently shorter oscillation period than the linear news. 
These results are of potential importance but, as already 
cautioned, at the current resolution it is not clear how 
accurate they are past u » 20. 



B. Azimuthal mode coupling 

We study azimuthal effects of mode coupling by pre- 
scribing initial data as an ingoing 2^22 pulse. Again 
we give the data in the form of Eq. (|3.3|1 with n = 3, 
%min = -56, and x max = .8, we vary the amplitude 
from A = 10~ 6 to A = .36 and we use a fixed grid of 
77 x 77 angular grid points and 181 radial grid points 
(larger grids are used when analyzing the cubic modes). 
The parity and reflection symmetries of the initial data 
arc preserved by the nonlinear evolution, so that the only 
non- vanishing amplitudes Ne m of the scattered radiation 
have even values of £ and even m > 0. However, in this 
case, both the © and ® polarization modes are excited. 
We decompose the resulting news function in terms of 
the spin-weighted functions 2^im up to t — 6. The non- 
vanishing quadratic modes present in the news were the 
2-R20, 2-R40, an d 2-R44 harmonics. The non- vanishing 
cubic modes were the 2-R42, 2-^62 and 2R66 harmonics. 

Fig. 1 121 shows the coefficient N22/A versus A. In linear 



FIG. 11: The news measured by an observer on the equator 
for A = 2.6 ■ lCT 1 and A = 2.6 • 10~ 6 . The amplitudes have 
been rescaled by a factor ~ 1/A to make comparison of the 
waveforms easier. 



theory the dependence of N22 on A is linear, and this 
is observed for the entire run when A < .1. Nonlinear 
amplification and phase shifts are only apparent for A — 
.36. 

Fig's. U21 El and graph the coefficients N 20 /A 2 , 
N40/A 2 , and N44/A 2 respectively. In each case there is 
a near perfect early time overlap between the curves at 
different amplitudes for A < .1. These nonlinear modes 
exhibit a quadratic dependence on A except for a late 
time amplification of the A = .36 curves in all modes and 
a late time amplification of the A = .1 curve of the N40 
mode. This late time amplification results from cubic and 
higher terms contributing to the modes. The very strong 
late-time behavior observed in N40 may result from the 
loss of resolution near r = 2M. 

Figures Hfjl 1171 and 1181 graph the coefficients N42/A 3 , 
NQ2/A 3 , and Nqq/A 3 . These modes show a cubic depen- 
dence on A, and were obtained by using a computational 
grid of 137 x 137 angular gridpoints and 331 radial grid- 
points (roughly twice the resolution of the previous runs). 
The amplitudes for these runs were A = .20, A = .26, and 
A = .36. 

The early time behavior of the {I = 2, m = 0) mode, 
unlike the behavior of the other nonlinear modes, is domi- 
nated by a low frequency component. However, hidden in 
this large signal is a higher frequency signal with roughly 
twice the frequency of the (£ = 2, m = 2) mode. Fig. 1191 
shows the Fourier decomposition for these modes. The 
{I = 2, m = 2) mode has a strong peak at u> w .33, the 
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FIG. 12: Azimuthal coupling. The rescaled coefficient N2 2 /A 
is plotted versus time. Note the near perfect linear depen- 
dence of N22 on A for A < .1. Nonlinear effects are only 
evident for the largest amplitude. 



(£ = 2, m = 0) mode has strong peaks atw» .80 and 
oj w .18, and both the (I = 4, m = 0) and the (t = 4, 
m = 4) modes have strong peaks at u> « .70. Given 
the large widths of the peaks, they are roughly consis- 
tent with the expected behavior of modes generated by 
quadratic response to an (£ — 2, m = 0) mode. On 
the other hand the frequency spectrum of the (£ = 4, 
m = 2) mode has a maximum atw« .63. The expected 
frequency of a cubic mode is either the frequency of the 
principal mode or three times the principal frequency. 
This drift toward larger than expected frequencies mir- 
rors the behavior of the (£ = 2, m = 0) mode. The 
(£ = 6, m = 2) mode has a maximum atw« 1.23 while 
the (I = 6, m = 6) mode has a maximum at u) ~ 1.02. 
These latter two results are roughly consistent with the 
expected behavior for cubic modes. 

Fig. [20| shows the Bondi news observed at q = p = .5 
(0 ~ 70 o ,(/) = 45°) for the A = .36, A = HT 1 , and 
v4 = 10~ 6 non-axisymmetric runs. In the linearized ap- 
proximation, the news in this angular direction is always 
imaginary, corresponding to pure © polarization. Non- 
linear effects couple the © and © modes. Initial data 
with A — .1 produce a significant © component with am- 
plitude roughly 28% of the © component; while initial 
data with A = .36 (at u = 20) produces an © component 
28% larger than the © component. Note that there is 
no significant nonlinear change in amplitude for the © 
component for A < .1. However, if a gravitational wave 
detector were not precisely oriented to measure the © 
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FIG. 13: Azimuthal coupling. The rescaled coefficient 
N20/A 2 is plotted versus time. Note the near perfect 
quadratic dependence of N20 on A for A < .1. Hidden within 
the low frequency signal at u < 15 there is a higher frequency 
mode as is apparent in Fig. 1191 The slight deviation of the 
A — 10~~ 6 curve is due to roundoff error. 



component, significant nonlinear amplification and phase 
shifting would arise from the superposition of the © and 
© components. In the axisymmetric case, all of the grav- 
itational radiation is in the © mode, and a gravitational 
wave detector would see the same nonlinear amplification 
and phase shifts regardless of orientation. In the present 
case the degree of nonlinear amplification and phase shift- 
ing depends on the orientation of the detector. 

We can gain insight into this problem by looking at the 
behavior of the real, spin-weight zero spherical harmon- 
ics Rim = QRim, (in a similar way as done in Sec. IV A"|l . 
For m even and positive, Re m is reflection symmetric 
about the q axis and has even parity with respect to the 
(q,p) coordinate system. For m even and negative, Rt m 
has even parity but is reflection antisymmetric about 
the q axis. Any product of positive, even m harmon- 
ics will contain these two symmetries and is therefore a 
sum of even, positive m, real spherical harmonics. Along 
with the reflection symmetry of Einstein's equations, this 
would seem to imply that data consisting of a positive, 
even m harmonic with even i would yield a Bondi news 
function consisting solely of positive, even m and even 
£ harmonics. In addition, we can predict the possible 
modes generated by an (£ — 2, m = 2) pulse and in what 
order they appear. For example, the quadratic modes 
produced by {£ = 2, m = 2) initial data are (£ = 2, 
m = 0), (£ = 4, m = 0), and (£ = 4, m = 4). Note 
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FIG. 14: Azimuthal coupling: The rescaled coefficient 
N40/A 2 is plotted versus time. Note the near perfect 
quadratic dependence of N40 on A at early times. The late- 
time deviation of the A = .36 curve may be due to the loss of 
resolution near r = 2M. The frequency is roughly twice that 
of N 2 2- 
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FIG. 15: Azimuthal coupling: The rescaled coefficient 
N44/A 2 is plotted versus time. Note the near perfect 
quadratic dependence of N44 on A for A < .1. The late 
time deviation is only significant for A > .1. The frequency 
is roughly twice that of -/V2 2 • 
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FIG. 16: Azimuthal coupling: The rescaled coefficient 
N42/A is plotted versus time. The overlap of the curves 
is reasonably good, but the frequency does not correspond to 
the expected behavior for a cubic mode. 
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FIG. 17: Azimuthal coupling: The rescaled coefficient 
NQ2/A 3 is plotted versus time. The overlap of the curves 
is reasonably good. Tests with higher resolution (269 x 269 
angular gridpoints) indicate that the poor early time scaling 
is due to truncation error. The frequency is roughly three 
times that of the -/V2 2 mode. 
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FIG. 18: Azimuthal coupling: The rescaled coefficient 
Nee/A 3 is plotted versus time. The overlap of the curves is 
reasonably good. The frequency is roughly three times that 
of the N22 mode. 

that there is no quadratic (£ = 2, to = 2) response. Thus 
quadratic terms in the Einstein equation do not affect 
the principal mode. Consequently there is no apparent 
nonlinear amplification or phase shift in this mode for 
A < .1 (see Fig. 11211 . Cubic terms can generate an (£ — 2, 
to = 2) mode so we expect those terms to produce non- 
linear amplification and phase shifts when they become 
significant. 

C. Generating catalogs of waveforms 

The results presented in See's. IVAl and IVBl suggest a 
method for efficiently producing a catalog of waveforms 
produced by a given profile of the initial data with am- 
plitude smaller than, say A = .1. First perform a run 
with the linearized code to generate the linear part of 
the waveform. Then perform a single nonlinear run at 
an amplitude of A = .1. One then subtracts off the lin- 
ear contribution to the various {£, to) modes and obtains 
the quadratic and cubic terms. The scaling with ampli- 
tude A of each of the (£, m) modes can be determined by 
the arguments given in See's. IVAl and IVBl and the net 
waveform for any A < .1 can be reconstructed. 

For example, consider the waveform generated from 
an axisymmctric {£ = 2, to — 0) input pulse with an 
amplitude A = .1. We first subtract off the linear part 
of the (£ = 2, m = 0) profile to get the quadratic part 
of the (£ = 2, m = 0) mode. We then assume a scaling 
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FIG. 19: Azimuthal coupling, The Fourier response (abso- 
lute value) to an A — .36 pulse is shown for the first few 
modes. The Fourier integrals were restricted to u < 15 and 
the resulting functions were rescaled to all have similar sizes. 
The dominant frequency of the N22 mode is uj « .33. The 
N2Q mode shows a strong low frequency peak at uj « .18 and 
a weaker peak at uj ~ .80. Interestingly, the transforms of 
-/V40 and A^4 4 are very similar with both having strong peaks 
at uj ~ .70. The -N42 mode (a cubic mode) shows an unex- 
pected strong peak at uj ~ .63. The 2 mode shows a strong 
peak at uj ~ 1.23 while the Nee mode has a strong peak at 
uj « 1.02. 

of A 2 for this part of the {£ = 2, m = 0) mode and the 
(I = 4, m = 0) mode, as well as a cubic scaling for the 
(I — 6, to = 0) mode. We can then produce the news for 
all A < .1 from a single nonlinear run. 

VI. CONCLUSION 

We have shown how a characteristic code can be em- 
ployed to investigate the nonlinear response of a black 
hole to the infall of gravitational wave energy. This prob- 
lem, which is of importance to the observation and inter- 
pretation of gravitational waves, can in this way be stud- 
ied with a mature, reliable evolution code. In this paper 
we have focused on the onset of nonlinear behavior in the 
waveform produced by the scattering of a pulse of radi- 
ation incident on a Schwarzschild black hole. However, 
given sufficient resolution, the same computational ap- 
proach could be extended to many other scenarios, such 
as the waveform emitted by the dispersion of a pulse of 
radiation propagating approximately along the r = 3M 
(unstable) orbit about a Schwarzschild black hole, as we 
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FIG. 20: Azimuthal ®, © coupling. The news is shown in 
the q — p = .5 observation direction and has been rescaled 
by 1/A. The rescaled <g> component of the yl = 1CP 1 and 
A — 10 -6 runs overlap. Nonlinear changes in amplitude and 
phase are visible in the <g) component of the A = .36 run. 

might expect from the gravitational perturbation associ- 
ated with a bounded distribution of matter in such an 
orbit. 

Besides the computation of accurate waveforms, our 
study reveals several features of qualitative importance: 

• I. The mode coupling amplitudes consistently scale 
as powers A n of the input amplitude A correspond- 
ing to the nonlinear order of the terms in the evo- 
lution equations which produce the mode. This al- 
lows much economy in producing a waveform cat- 
alog. Given the order n associated with a given 
mode generation, the response to any input am- 
plitude A can be obtained from the response to a 
single reference amplitude Aa. 

• II. The frequency response has similar behavior but 
in a less consistent way. The dominant frequencies 
produced by mode coupling are in the approximate 
range of the quasinormal frequency of the input 
mode and the expected sums and difference fre- 
quencies generated by the order of nonlincarity. 

• III Large phase shifts, ranging up 15% in a half 
cycle relative to the linearized waveform, are ex- 
hibited in the news function obtained by the super- 
position of all output modes, i.e. in the waveform 
of observational significance. These phase shifts, 
which are important for design of signal extraction 



templates, arise in an erratic way from superposing 
modes with different oscillation frequencies. This 
furnishes another strong argument for going be- 
yond the linearized approximation in designing a 
waveform catalog for signal extraction. 

• IV Besides the nonlinear generation of harmonic 
modes absent in the initial data, there is also 
a stronger than linear generation of gravitational 
wave output. This provides a potential mechanism 
for enhancing the strength of the gravitational ra- 
diation produced during, say, the merger phase of 
a binary inspiral above the strength predicted in 
linearized theory. 

• V In the non-axisymmetric case, there is also con- 
siderable generation of radiation in polarization 
states not present in the linearized approximation. 
In our simulations, input amplitudes in the range 
A = .1 to A = .36 lead to nonlinear generation of a 
© component which is of the same order of magni- 
tude as the <g> component (which would be the sole 
component according to linearized theory). As a re- 
sult, significant nonlinear amplification and phase 
shifting of the waveform can be observed depending 
on the orientation of a gravitational wave detector. 

As already noted by Papadopoulos in his work on ax- 
isymmetric mode coupling |2lj| . these effects arise from 
three types of nonlinearity: (i) Modification of the light 
cone structure governing the principal part of the equa- 
tions and hence the propagation of signals; (ii) Modula- 
tion of the Schwarzschild potential by the introduction of 
an angular dependent "mass aspect" ; and (iii) Quadratic 
and higher order terms in the evolution equations which 
couple the modes. 

Although Papadopoulos studied nonlinear mode gen- 
eration produced by an outgoing pulse, as opposed to the 
case of an ingoing pulse studied here, these same factors 
are in play and it is not surprising that both studies have 
common features. In both cases, the major nonlinear ef- 
fects arise in the region near r = 3M. Analogs of items I 
- IV above are all apparent in Papadopoulos 's results. At 
the finite difference level, both codes respect the reflec- 
tion symmetry inherent in Einstein's equations and ex- 
hibit the corresponding selection rules arising from parity 
considerations. In the axisymmetric case considered by 
Papadopoulos, this forbids the nonlinear generation of a 
© mode from a <g> mode, as in item V above. 

The evolution along ingoing null hypersurfaces in the 
work of Papadopoulos and the evolution along outgoing 
null hypersurfaces in the present work have complemen- 
tary numerical features. The grid based upon ingoing 
null hypersurfaces avoids the difficulty depicted in Fig. 
in resolving effects close to r = 2M with a grid based 
upon outgoing null hypersurfaces. The outgoing code 
would require some form of mesh refinement in order 
to resolve the quasinormal ringdown for as many cycles 
as achieved by Papadopoulos. However, the outgoing 
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code avoids the late time caustic formation noted in Pa- 
padopoulos's work, as well as the gauge ambiguity and 
backscattering which complicate the extraction of wave- 
forms on a finite worldtube. An attractive option is to 
combine the best features of both codes by matching an 
interior evolution based upon ingoing null hypersurfaces 
to an exterior evolution based upon outgoing null hyper- 
surfaces, as implemented in 36] for spherically symmetric 
Einstein-Klein-Gordon waves. 

The waveform of relevance to gravitational wave as- 
tronomy is the superposition of modes with different fre- 
quency compositions and angular dependence. Although 
this waveform results from a complicated nonlinear pro- 
cessing of the input signal, which varies with choice of 
observation angle, we have shown that the response of 
the individual modes to an input signal of arbitrary am- 
plitude can be obtained by scaling the response to an 
input of standard reference amplitude. This offers an 
economical approach to preparing a waveform catalog. 

Work is in progress pursuing the above projects, as 
well as extending the treatment to spinning black holes. 
In concert with other approaches to compute nonlinear 
waveforms, it is hoped that robust features of the gravi- 
tational waves produced by highly distorted black holes 
will be discovered which could be exploited in data anal- 
ysis efforts. 
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APPENDIX A: SPIN WEIGHTED SPHERICAL 
HARMONICS 

The complex stereographic coordinate z = q + ip cov- 
ers the sphere with two patches which overlap in a region 
containing the equator. In the north patch, the stereo- 
graphic coordinate is related to standard spherical coor- 
dinates (6, 4>) by zn — tan |e Ic ^; and in the south patch, 



by zs = cot | e~^. In the overlap region, the coordinates 
are related by z$ = 1/zjv- 

In the nor th p atch, the ordinary spherical harmonics 
are given by |38| 



min(e,l+m) 
n— max(0,m) 



n\{£ + m — n)\{£ — n)\(n — m)\ 



pi 



VT+21 „, (£-m)\(£ + m)\ 
x — — — 1\\ 



4tt 



(Al) 



where P = 1 + q 2 + p 2 . We give the explicit formulae 
and conventions for the spin- weighted spherical harmon- 
ics used in our calculations in order to avoid confusion 
with various conventions found in the literature. The 
spin- weighted spherical harmonics arc defined by 



(e+s)\ 



d s Y tm , s>0, (A2) 



s Y im = (-l)y^B- s Y lm , s<0 (A3) 

in terms of the spin-weight raising and lowering operators 
3 and 3 which act on a spin- weight s function / according 
to 

3/ = P 1 -'d s (fP a ), Bf = P 1+s d z (fp- s ). (A4) 
They are given in the north patch by 

min(£-t- s,£+m) 



v V ( iv VTTYe 

^ {I + znznY 

n-max(0,s-fm) 



-N^N 



n\{£ + m — n)\(£ + s — n)!(n — m — s)\ 

x l (£-m)\(£ + m)\(e-s)\(e + s y^ (A5) 

The value fs(zs) of any spin- weight s function / in the 
south stereographic coordinates is given in terms of its 
value Jn{zn) in north stereographic coordinates by [3(| 



fs(zs) 



ZS_ 

zs 



fif(l/zs). 



(A6) 



Consequently, the spin-weighted spherical harmonics in 
the south stereographic patch are given by 



min(£+s,£-(-m) 



,Y S 



s o vm 



V ( L) s +" v 1 + M 
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The spin-weighted spherical harmonics obey the or- 
thogonality condition 



(A8) 
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where the integration is over the entire solid angle £1 = An 
of the unit sphere. The integral l|A8[) can be re-expressed 
as 



s^£m sYg' m' dQ — 



+ 



4( s Yn lm){ s Yn I'm') 

n (± + q 2 N +p 2 N ) 2 

HsYsi m )( s Y s 

V nv ) 



dqN dp N 



(l + ql+p 2 s ) 



2 \2 



dq s dps (A9) 



where the subscripts N and S refer to the north and south 
patches respectively and J N and J s denote integration 
over the corresponding hemispheres. 

Rather than decompose J and the news in terms of 
the sYlm we decompose these functions in terms of the 
sRim which are defined by 

sRlm = -j= [ s Yl m + {-l) m s Yl- m } for TO > 



sRim - j= 

sRia = sYgo 



[{-l) m s Y em - s Y e _ m ] forTO<0 

(A10) 



These superpositions of the s Yi m obey the orthogonality 
condition 



,Rim sRe'm' dCl = 5ee'S n 



(All) 



The advantage of using the s Rg m as basis functions is 
that in the linear regime there is no coupling between 
different to modes. The news calculation contains terms 
of the form 9 2 (3 2 J + 3 2 J) which introduces linear order 
mode coupling in the s Yg m basis on account of 

S 2 2^m + B 2 2 Y em ~ X(Y tm ) ~ Y lm + {-l) rn Y e - m . 

However, since $l(Ri m ) = Re m there is no spurious m 
mode coupling in the s Rim basis. 

Any smooth spin-weight s function F can be decom- 
posed into the corresponding spin- weighted spherical har- 
monics according to 



F 



Ftrn sR 



t ra i 



where 



Fim — <t> F s Rf m dQ. 



(A12) 



(A13) 



The decomposition in Eq. (1 A 1211 defines a spin-zero po- 
tential / for the spin-weighted function F given by 



£m 

/ = -i s £ 



Cm 



'{£-s)\ 
(£ + s)\ 

l (e + s)\ 

(£-s)\ 



FlmRlrn for S > 0, (A14) 



FimRim for s < 0, (A15) 



where Rg m = oRi m , F — 9 s / for s > 0, and F = 5 s f 
for s < 0. Of particular relevance are the coefficients of 



the various modes of the Bondi news function given by 



T N ff 2 Rn £m dqN dpN 



Ntm = j N 2 R em dQ 

= [ 1 

Jn + 

4 
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APPENDIX B: MODIFICATIONS TO THE NEWS 
MODULE 

In order to calculate the news function in an inertial 
Bondi gauge one needs to track the phase angle which ro- 
tates the complex dyad defined using the coordinates of 
the PITT null code into the complex unit sphere dyad de- 
fined with respect to an inertial Bondi coordinate system 
on 1+ (see Eq's (32)- (36) of [3). This phase, which 
we denote by e , obeys a hyperbolic partial differen- 
tial equation in the non-inertial coordinates of the code 
(Eq. (36) of 0) 



(d u + L A d A )6 



1 / J iU J J (UBJ + UdJ) 
2^\K + 1 + 2(^ + 1) 



JW + KW + 2Uz 



(Bl) 



The solution of Eci . IBTI rcci uires consistent boundary data 
at the edges of the stereographic patches. Because 5 is 
not a pure spin-weighted function (as is evident in the Uz 
term in Eq. IjBll) - ). the necessary cross-patch interpolation 
rules are complicated. However, we have found that the 
computation of S can be simplified by recasting Ea. lBll as 
an ordinary differential equation along the characteristics 
of the equation, which are the null generators of 2 + . In 
the inertial Bondi coordinates, this ordinary differential 
equation takes the form 



dS 
du 



J U J J(UBJ+Udl) 



2 \K + 1 
JW + KW 



2(K + 1) 
2Uz) 



(B2) 



where z and u are the non-inertial coordinates. In this 
formulation no boundary data are required. Eq. <|B2|) is 
integrated to second order accuracy via 



rn+l 



1 



du (RHS™ + RHS 



(B3) 



where RHS is the right-hand-side of Eq. i|B2|l , the index 
i labels the null generators of T + , and the index n indi- 
cates the time level. (Although the above modification 
has been implemented in the present code, it is not neces- 
sary when considering axisymmetric spacetimes. In that 
case Q(Uz) = and S behave as ordinary spin-weight 
functions.) For more details on the numerical implemen- 
tation, see 391. 
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